R Notebook

library(readxl)

setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
feature_id <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "feature_id")
features <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "features")
metadata <- read_excel("metabolomics_microbiome_fc_2023.xlsx", sheet = "metadata")

library(dplyr)
library(tibble)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")

diff_all <- features %>% column_to_rownames("features") %>% t() %>% as.data.frame() %>% rownames_to_column(var="sample_ID") %>% merge(metadata, by="sample_ID")
matrix_diff <- diff_all %>% remove_rownames() %>% column_to_rownames(var="sample_ID")
matrix_diff=as.matrix(matrix_diff)
class(matrix_diff) <- "numeric"

Warning: NAs introduced by coercion

matrix_diff[is.na(matrix_diff)] = 0

#ANNOTATION COLUMNS
anno_row <- diff_all %>% 
  dplyr::select(Treatment,Time_tx)

rownames(anno_row)= rownames(as.data.frame(matrix_diff))

anno_row <- as.data.frame(anno_row)

anno_color <- list(Treatment = c(Antibiotic = "#DF8F44FF", Control = "#374E55FF"),
                   Time_tx = c(`Week 1` = "#EE0000FF", `Day -1` = "#3B4992FF", `Week 9` = "#008B45FF"))

library(pheatmap)
library(viridis) #color pallet, it's optional
library(RColorBrewer)

#HEATMAP RA USING ln color pallet RBrewer
heatmap_metIDs <- pheatmap(
    #mat               = t(matrix_diff),
    mat               = log10(t(matrix_diff+0.0000001)),
  border_color      = NA,
  show_colnames     = T,
  show_rownames     = T,
  angle_col = 90,
  drop_levels       = TRUE,
  fontsize_col = 4,
  fontsize_row = 5,
  fontsize          = 14,
  color             = brewer.pal(9,"Reds"),
  number_color = NA,
 annotation_col  = anno_row,
  annotation_colors = anno_color,
  annotation_names_col = T,
  annotation_names_row = T,
  cluster_cols = T,
  cluster_rows = T,
 clustering_method = "ward.D",
  gaps_row = FALSE
 )

Warning: NaNs produced

heatmap_metIDs

setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
matrix_features <- read_excel("metabolomics_microbiome_fc.xlsx", sheet = "features")%>% remove_rownames %>% column_to_rownames(var="features")
matrix_features=as.matrix(matrix_features)
class(matrix_diff) <- "numeric"
matrix_diff[is.na(matrix_diff)] = 0
hc <- hclust(dist(matrix_features), "ward.D")
plot(hc, cex =.5)

hc_groups <- as.data.frame(cutree(hc, k = 30)) %>% tibble::rownames_to_column("features") %>% rename(hc_group = 2)

hc_groups %>% count(hc_group)

#Number of clusters

library("NbClust")
nb <- NbClust(diss = dist(matrix_features), distance = NULL, min.nc = 2,
        max.nc = 50, method = "ward.D",index = "silhouette")

 Only frey, mcclain, cindex, sihouette and dunn can be computed. To compute the other indices, data matrix is needed 

library(janitor)
library(tidyr)
diff <- as.data.frame(matrix_features) %>% rownames_to_column("features") 
df_hc_merge <- merge(diff, hc_groups, by = "features")
df_hc <- df_hc_merge %>% gather(sample_ID, abundance, -c(features,hc_group)) 
#metadata <- metadata %>% rename(sample_ID = ID)
df_hc_meta <- merge(df_hc, metadata, by="sample_ID")
df_hc_meta$abundance <- as.numeric(as.character(df_hc_meta$abundance))
df_hc_meta_type <- merge(df_hc_meta, feature_id, by="features")

hc2 <- df_hc_meta_type %>% 
  filter(hc_group==2, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,15)
hc2

hc4 <- df_hc_meta_type %>% 
  filter(hc_group==4, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,5)
hc4

hc5 <- df_hc_meta_type %>% 
  filter(hc_group==5, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc5

hc6 <- df_hc_meta_type %>% 
  filter(hc_group==6, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,5)
hc6

hc7 <- df_hc_meta_type %>% 
  filter(hc_group==7, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc7

hc8 <- df_hc_meta_type %>% 
  filter(hc_group==8, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc8

hc9 <- df_hc_meta_type %>% 
  filter(hc_group==9, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,10)
hc9

hc10 <- df_hc_meta_type %>% 
  filter(hc_group==10, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc10

hc11 <- df_hc_meta_type %>% 
  filter(hc_group==11, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ 
  theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc11

hc12 <- df_hc_meta_type %>% 
  filter(hc_group==12, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ 
  theme(legend.text=element_text(size=7))  + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc12

hc14 <- df_hc_meta_type %>% 
  filter(hc_group==14, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx",panel.labs = list(Time_tx=c("Lactation","Dry-off","Fresh cows"))) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+  guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=c("#DF8F44FF","#B24745FF")) + scale_color_manual(values=c("#DF8F44FF","#B24745FF"))
hc14

hc15 <- df_hc_meta_type %>% 
  filter(hc_group==15, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ 
  theme(legend.text=element_text(size=7))  + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ylim(0,10)
hc15

hc28 <- df_hc_meta_type %>% 
  filter(hc_group==14, Time_tx%in%c("Day -1","Week 1","Week 9")) %>% 
  ggboxplot(x="features",y="abundance",fill = "Type",color = "Type", ylab = "Mean proportion", alpha = 0.5, notch = F,legend = "right",facet.by = "Time_tx", palette = get_palette(palette = "jama",5)) + theme(axis.ticks.x = element_blank(), axis.title.x = element_blank())+ theme(legend.text=element_text(size=5)) + guides(fill=guide_legend(ncol=1),color=guide_legend(ncol=1))+ 
  theme(legend.text=element_text(size=7)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
hc28

setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")
#ggsave(plot=hc1,"hc1_multiomics.png",width = 16, height = 8)
ggsave(plot=hc2,"hc2_multiomics.png",width = 30, height = 20)
#ggsave(plot=hc3,"hc3_multiomics.png",width = 20, height = 20)
#ggsave(plot=hc4,"hc4_multiomics.png",width = 16, height = 8)
ggsave(plot=hc5,"hc5_multiomics.png",width = 30, height = 20)
#ggsave(plot=hc14,"hc14_multiomics.png",width = 16, height = 20)
#ggsave(plot=hc15,"hc15_multiomics.png",width = 16, height = 20)

setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")

hc_bars <- ggarrange(hc1,hc14,hc3,labels = c("A","B","C"),nrow = 3, ncol = 1)
ggsave(plot=hc_bars,"hc_bars_fc_multiomics.png",width = 15, height = 15)